####################################################
####################################################
###### Code for Hannah, A. Lee and Daniel J. Mallinson. Forthcoming. "Defiant Innovation: 
###### The Adoption of Medical Marijuana Laws 
###### in the U.S. States" Policy Studies Journal. 
###### 27 March 2017
###### lee.hannah@wright.edu 
####################################################
####################################################


rm(list=ls()) #clear workspace

#Read adoption model data 
library(foreign)
# setwd()
raw.data <- read.dta("MM_Duration_Final_Revision.dta") #Raw data
data <- raw.data
head(data)

data <- data[which(data$year > 1995),] #Drop 1995 observations 


####################################################
###### 1. CODE TO PRODUCE RELATIVE IDEOLOGY  #######
####################################################

df <- as.data.frame(matrix(NA, nrow=1, ncol=ncol(data)+1)) #Create blank dataset for appending
names(df) <- c(colnames(data), "ideology_relative_berry")

#for(i in 1:nrow(unique.policies)){
  #use.policy <- unique.policies[i,]
  #use.data <- data[which(data$policy==use.policy),]
  vars <- c("state", "year","fip", "citi6013", "medical_marijuana")
  adopts <- data[which(data$medical_marijuana==1),]
  adopts <- adopts[vars]
  for(j in 1:nrow(data)){
    use.obs <- data[j,]
    if(use.obs$year<1997){
      ideo.first <- data[which(data$year==min(data$year)),]
      ideo <- mean(ideo.first$citi6013)
      use.obs$ideology_relative_berry <- abs(use.obs$citi6013 - ideo)
      df <- rbind(df, use.obs)}
    else{
      ideo.adopts <- adopts[which(adopts$year<use.obs$year),]
      use.obs$ideology_relative_berry <- (sum(abs(use.obs$citi6013 - ideo.adopts$citi6013))/nrow(ideo.adopts))
      df <- rbind(df, use.obs)}
  }
#}

head(df)
df <- df[-1,] #Delete row of null values

head(df)

data <- df

write.csv(df, file="data_with_ideo_6.csv", row.names=FALSE)





####################################################
###### 2. REPRODUCTION CODE FOR FINDINGS     #######
#######            IN PAPER                  #######
####################################################


setwd("PSJ Final") ## Set directory to find file

library(stats)
library(foreign)

rm(list=ls()) #clear workspace

raw.data <- read.csv("data_with_ideo_6.csv")

#Pull only variables used in the analysis and eliminate missingness. This is necessary for the calculation of clustered standard errors if some states are dropped because of missingness. 
variables <- c("fip", "state", "year", "medical_marijuana", "evan_rate", "init_avail", "bush_admin", "obama_admin","prop_neighbor", "ideology_relative_berry","cancer_rt_sl", "glaucoma_rt", "lp_avg", "marijuana_kgs", "marijuana_ted", "fiscal_health")


#data <- na.omit(raw.data[variables])
data <- raw.data

colnames(data)
head(data)
#fix(data)

## Recscale marijuana into 10,000 kgs
data$marijuana_kgs_10k = (data$marijuana_kgs/10000)
data$marijuana_lag_kgs_10k = (data$marijuana_lag/10000)

#Rescale glaucoma rate into percentage
data$glaucoma_pct = (data$glaucoma_rt*100)
data$time_adopt = (data$year-1996)

###################################################
###### CODE TO PRODUCE TABLE 2 COEFFICIENTS  ######
###################################################


model.full <- glm(medical_marijuana ~ marijuana_kgs_10k + bush_admin + obama_admin + init_avail + lp_avg + evan_rate + citi6013 + fiscal_health + marijuana_ted + glaucoma_pct + cancer_rt_sl + prop_neighbor + ideology_relative_berry + time_adopt, data=data, family=binomial)
summary(model.full)

###################################################
########   TABLE 2 Model Fit Statistics   ######### 
###################################################
library(aod)
wald.test(b=coef(model.full), Sigma=vcov(model.full), Terms=2:15)
logLik(model.full)


###################################################
###########    TABLE 2 Robust SEs       ########### 
###################################################

library(sandwich)
library(zoo)
use.model <- model.full

fc <- data$fip
m <- length(unique(fc))
k <- length(coef(use.model))
vcov <- vcov(use.model)

u <- estfun(use.model)
u.clust <- matrix(NA, nrow=m, ncol=k)
for(j in 1:k){
	u.clust[,j] <- tapply(u[,j], fc, sum)
}

cl.vcov <- vcov %*% ((m/(m-1)) * t(u.clust) %*% (u.clust)) %*% + vcov

library(lmtest)
coeftest(use.model, cl.vcov)

###################################################
###########    TABLE 2 EFFECTS SIZE     ########### 
###################################################
# Rerun model.full

library(lmtest)
coeftest(use.model, cl.vcov)

#Marginal Effects
library(mfx)
library(betareg)
logitmfx(medical_marijuana ~ marijuana_kgs_10k + bush_admin + obama_admin + init_avail + lp_avg + evan_rate + citi6013 +fiscal_health + marijuana_ted + glaucoma_pct + cancer_rt_sl + prop_neighbor + ideology_relative_berry + time_adopt, data=data, robust=TRUE, atmean=FALSE)

#Odds Ratios
logitor(medical_marijuana ~ marijuana_kgs_10k + bush_admin + obama_admin + init_avail + lp_avg + evan_rate + citi6013 + fiscal_health + marijuana_ted + glaucoma_pct + cancer_rt_sl + prop_neighbor + ideology_relative_berry + time_adopt, data=data, robust=TRUE)

X <- summary(use.model)$coef

or <- as.data.frame(X[,1:2])
or$odds <- exp(or[,1])
or$odds.lowerci <- exp(or[,1] - (1.96*or[,2]))
or$odds.upperci <- exp(or[,1] + (1.96*or[,2]))
or

###################################################
######   CALCULATE % CHANGE IN ODDS (TEXT)   ###### 
###################################################

#Percent change in odds
#Bush
100*(0.021-1) #-97.9

#Obama
100*(0.008-1) #-99.2

#Initiative
100*(3.03-1) #203

# Ideology
100*(1.097 - 1) #9.7
sd(data$citi6013) #14.99

9.7*14.99 #145.403

# Evangelical
100*(0.989-1) #-1.1
sd(data$evan_rate) #117.7849
-1.1*117.7849 #-129.56

#Relative Ideology
100*(0.881-1) #-11.9
sd(data$ideology_relative_berry) #8.171935 
-11.9*8.171935 #-97.25






